function [t,R,V,R_maneuver,t_maneuver,delVs,...
    maneuverstart_utc,maneuverfinish_utc,maneuvertype,...
    ManeuverMass,ManeuverDuration]...
    = SK_GEO_LT_tolerance_simulation...
    (sim_day,tol_lat,tol_long,...
    target_long,station_keeping_duration_max_hr,mass,isp,T)
%%
global const Cnm Snm AuxParam eopdata swdata SOLdata DTCdata APdata PC

SAT_Const
constants
load DE436Coeff.mat
PC = DE436Coeff;

% read Earth gravity field coefficients
Cnm = zeros(181,181);
Snm = zeros(181,181);
fid = fopen('GGM03S.txt','r');
for n=0:180
    for m=0:n
        temp = fscanf(fid,'%d %d %f %f %f %f',[6 1]);
        Cnm(n+1,m+1) = temp(3);
        Snm(n+1,m+1) = temp(4);
    end
end
fclose(fid);

% read Earth orientation parameters
fid = fopen('eop19620101.txt','r');
%  ----------------------------------------------------------------------------------------------------
% |  Date    MJD      x         y       UT1-UTC      LOD       dPsi    dEpsilon     dX        dY    DAT
% |(0h UTC)           "         "          s          s          "        "          "         "     s
%  ----------------------------------------------------------------------------------------------------
eopdata = fscanf(fid,'%i %d %d %i %f %f %f %f %f %f %f %f %i',[13 inf]);
fclose(fid);

% read space weather data
fid = fopen('sw19571001.txt','r');
%  ---------------------------------------------------------------------------------------------------------------------------------
% |                                                                                             Adj     Adj   Adj   Obs   Obs   Obs
% | yy mm dd BSRN ND Kp Kp Kp Kp Kp Kp Kp Kp Sum Ap  Ap  Ap  Ap  Ap  Ap  Ap  Ap  Avg Cp C9 ISN F10.7 Q Ctr81 Lst81 F10.7 Ctr81 Lst81
%  ---------------------------------------------------------------------------------------------------------------------------------
swdata = fscanf(fid,'%4i %3d %3d %5i %3i %3i %3i %3i %3i %3i %3i %3i %3i %4i %4i %4i %4i %4i %4i %4i %4i %4i %4i %4f %2i %4i %6f %2i %6f %6f %6f %6f %6f',[33 inf]);
fclose(fid);

% read solar storm indices
fid = fopen('SOLFSMY.txt','r');
%  ------------------------------------------------------------------------
% | YYYY DDD   JulianDay  F10   F81c  S10   S81c  M10   M81c  Y10   Y81c
%  ------------------------------------------------------------------------
SOLdata = fscanf(fid,'%d %d %f %f %f %f %f %f %f %f %f',[11 inf]);
fclose(fid);

% read geomagnetic storm indices
fid = fopen('DTCFILE.txt','r');
%  ------------------------------------------------------------------------
% | DTC YYYY DDD   DTC1 to DTC24
%  ------------------------------------------------------------------------
DTCdata = fscanf(fid,'%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d',[26 inf]);
fclose(fid);

% read geomagnetic storm indices
fid = fopen('SOLRESAP.txt','r');
%  ------------------------------------------------------------------------
% | YYYY DDD  F10 F10B Ap1 to Ap8
%  ------------------------------------------------------------------------
APdata = fscanf(fid,'%d %d %f %f %f %f %f %f %f %f %f',[12 inf]);
fclose(fid);

% model parameters
AuxParam = struct('Mjd_UTC',0,'area_solar',0,'area_drag',0,'mass',0,'Cr',0,...
    'Cd',0,'n',0,'m',0,'sun',0,'moon',0,'sRad',0,'drag',0,...
    'planets',0,'SolidEarthTides',0,'OceanTides',0,'Relativity',0);
%% Do not Touch
% epoch state (Envisat)
fid = fopen('InitialCondition_data.txt','r');
tline = fgetl(fid);
year = str2num(tline(1:4));
mon = str2num(tline(6:7));
day = str2num(tline(9:10));
hour = str2num(tline(12:13));
min = str2num(tline(15:16));
sec = str2num(tline(18:19));
UTC = [year,mon,day,hour,min,sec];
Y0 = zeros(1,6);
for ii=1:6
    tline = fgetl(fid);
    Y0(ii) = str2num(tline);
end
tline = fgetl(fid);
AuxParam.area_solar = str2num(tline(49:end));
tline = fgetl(fid);
AuxParam.area_drag = str2num(tline(38:end));
tline = fgetl(fid);
AuxParam.mass = str2num(tline(19:end));
tline = fgetl(fid);
AuxParam.Cr = str2num(tline(5:end));
tline = fgetl(fid);
AuxParam.Cd = str2num(tline(5:end));
fclose(fid);

% epoch
Mjd_UTC = Mjday(year, mon, day, hour, min, sec);

AuxParam.Mjd_UTC = Mjd_UTC;
AuxParam.n       = 70;
AuxParam.m       = 70;
AuxParam.sun     = 1;
AuxParam.moon    = 1;
AuxParam.planets = 0;
AuxParam.sRad    = 0;
AuxParam.drag    = 0;
AuxParam.SolidEarthTides = 0;
AuxParam.OceanTides = 0;
AuxParam.Relativity = 0;
Mjd0   = Mjd_UTC;
Step   = 60;   % [s]
We = (2*pi / 86164.1) * [0;0;1];
%[] rotation speed of Earth

%% Change from here


simulation_day = sim_day;
%[] Total simulation days

target_lat = 0;
%[] Longitude and latitude tolerance.

d_day = 5;
%[] How many days of integration per while loop ?

N_Step = d_day*24*60*60/Step;
%[] Change d_day to step/

R_target_ecef = Eci2Ecef(UTC)*Y0(1:3)';
V_target_ecef = Eci2Ecef(UTC)*Y0(4:6)' - cross(We,R_target_ecef);
S_target_ecef = [R_target_ecef;V_target_ecef]/1000;
%[]State of the perfect location in ECEF frame.

MAX_STATIONKEEPING_DURATION_HOUR = station_keeping_duration_max_hr;
%[] Search space of station keeping duration

t_save = [];
R_save = [];
V_save = [];
R_maneuver_save = [];
t_maneuver_save = [];
Number_of_maneuver = 0;
maneuverstart_utc = [];
maneuverfinish_utc = [];
maneuvertype = {};
delVs = [];
ManeuverMass = [];
ManeuverDuration = [];

cexh = isp*9.81;
%[] Create  read/write-able memory

currentTime = 0;
% Current Time 

t0 = 0;
%Initial time from given UTC

target_accel = T/mass;



while currentTime < simulation_day*24*60*60
    
    % Propagation
    Eph = Ephemeris(Y0, N_Step, Step,t0);
    Eph = Eph';
    t = Eph(1,:);
    R = Eph(2:4,:)/1000;
    V = Eph(5:7,:)/1000;
    n = length(t);
    [latitude,longitude] = ECI_Position2LatLong(t,R,UTC);
    binary_latitude_check = or((latitude > target_lat+tol_lat),(latitude < target_lat-tol_lat));
    binary_longitude_check = or((longitude > target_long+tol_long),(longitude < target_long-tol_long));
    
    
    % Check if outside of tolerance
    check1 = strfind(binary_latitude_check,[0,1]);
    if length(check1)~=1
        if isempty(check1)==1
            check1 = inf;
        end
        check1 = check1(1);
    end
    check2 = strfind(binary_longitude_check,[0,1]);
    if length(check2)~=1
        if isempty(check2) == 1
            check2 = inf;
        end
        check2 = check2(1);
    end
    
    
    
    % Skip
    if and(check1 == inf,check2 == inf)
        t_save = [t_save,t];
        R_save = [R_save,R];
        V_save = [V_save,V];
        Y0 = [R_save(:,end)*1000;V_save(:,end)*1000];
        t0 = t_save(end);
        currentTime = t0(end);
        continue
    
        
    % Perform maneuver for latitude tolerance
    elseif check1(1) < check2(1)
        Mjd_UTC_new = Mjday(year, mon, day, hour, min, sec+t(check1(1)));
        UTC_new = [year, mon, day, hour, min, sec+t(check1(1))];
        R_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(1:3);
        V_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(4:6) + cross(We,R_target_eci);
        S_target_eci = [R_target_eci;V_target_eci];
        t_crossing = t(check1(1));
        [R_maneuver,V_maneuver,t_maneuver,delV] = Run_StationKeeping_LT([R(:,check1(1));V(:,check1(1))],S_target_eci,S_target_ecef,...
                                                        t_crossing,UTC,target_accel,mass,isp,T);
        t_save = [t_save,t(:,2:check1(1)),t_maneuver];
        R_save = [R_save,R(:,2:check1(1)),R_maneuver];
        V_save = [V_save,V(:,2:check1(1)),V_maneuver];
        R_maneuver_save = [R_maneuver_save,R_maneuver];
        t_maneuver_save = [t_maneuver_save,t_maneuver];
        Y0 = [R_save(:,end);V_save(:,end)]*1000;
        t0 = t_save(end);
        currentTime = t0(end);
        Number_of_maneuver = Number_of_maneuver+1;
        maneuverstart_utc = [maneuverstart_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(1)]))];
        maneuverfinish_utc = [maneuverfinish_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(end)]))];
        maneuvertype = [maneuvertype;'North-South'];
        delVs = [delVs;delV];
        ManeuverMass = [ManeuverMass;mass*exp(delV/cexh)-mass];
        ManeuverDuration = [ManeuverDuration;etime(maneuverfinish_utc(end,:),maneuverstart_utc(end,:))];
        continue
        
        
    % Perform maneuver for longitude tolerance
    else
        Mjd_UTC_new = Mjday(year, mon, day, hour, min, sec+t(check2(1)));
        UTC_new = [year, mon, day, hour, min, sec+t(check2(1))];
        R_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(1:3);
        V_target_eci=transpose(Eci2Ecef(UTC_new))*S_target_ecef(4:6) + cross(We,R_target_eci);
        S_target_eci = [R_target_eci;V_target_eci];
        t_crossing = t(check2(1));
        [R_maneuver,V_maneuver,t_maneuver,delV] = Run_StationKeeping_LT([R(:,check2(1));V(:,check2(1))],S_target_eci,S_target_ecef,...
                                                        t_crossing,UTC,target_accel,mass,isp,T);
        t_save = [t_save,t(:,2:check2(1)),t_maneuver];
        R_save = [R_save,R(:,2:check2(1)),R_maneuver];
        V_save = [V_save,V(:,2:check2(1)),V_maneuver];
        R_maneuver_save = [R_maneuver_save,R_maneuver];
        t_maneuver_save = [t_maneuver_save,t_maneuver];
        Y0 = [R_save(:,end);V_save(:,end)]*1000;
        t0 = t_save(end);
        currentTime = t0(end);
        Number_of_maneuver = Number_of_maneuver+1;
        maneuverstart_utc = [maneuverstart_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(1)]))];
        maneuverfinish_utc = [maneuverfinish_utc;CalendarDate(JulianDate([UTC(1:5),UTC(6)+t_maneuver(end)]))];
        maneuvertype = [maneuvertype;'East-West'];
        delVs = [delVs;delV];
        ManeuverMass = [ManeuverMass;mass*exp(delV/cexh)-mass];
        ManeuverDuration = [ManeuverDuration;etime(maneuverfinish_utc(end,:),maneuverstart_utc(end,:))];

        continue
    end
end

t = t_save;
R = R_save;
V = V_save;
R_maneuver = R_maneuver_save;
t_maneuver = t_maneuver_save;

end

